function B = PMSqu3D(M,ao,l,position)
%func:@使用解析模型计算条形磁铁的磁感应强度(B)
%inpt:@M:剩磁场
%     @(ao,l):磁铁参数
%     @(xd,yd,zd)被计算的位置
%oupt:@B = [Bx;By;Bz]
%note:底面可以更改为非正方形;
%date:2019/12/18 by xiyin,li in ChongQing University;
%@Copyright 2019 Xiyin,li.
xd = position(1);
yd = position(2);
zd = position(3);
nu0 = 4*pi*10^(-7);
ro = ao/2;
Bx1 = dblquad(@(x,y)PMSquFunc(x,y,l/2,xd,yd,zd,1),-ro,ro,-ro,ro,1e-4)/(4*pi)*M;
Bx2 = dblquad(@(x,y)PMSquFunc(x,y,-l/2,xd,yd,zd,1),-ro,ro,-ro,ro,1e-4)/(4*pi)*M;
By1 = dblquad(@(x,y)PMSquFunc(x,y,l/2,xd,yd,zd,2),-ro,ro,-ro,ro,1e-4)/(4*pi)*M;
By2 = dblquad(@(x,y)PMSquFunc(x,y,-l/2,xd,yd,zd,2),-ro,ro,-ro,ro,1e-4)/(4*pi)*M;
Bz1 = dblquad(@(x,y)PMSquFunc(x,y,l/2,xd,yd,zd,3),-ro,ro,-ro,ro,1e-4)/(4*pi)*M;
Bz2 = dblquad(@(x,y)PMSquFunc(x,y,-l/2,xd,yd,zd,3),-ro,ro,-ro,ro,1e-4)/(4*pi)*M;
B = [Bx1-Bx2;By1-By2;Bz1-Bz2];